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Abstract 

To reduce scanning time and/or improve spatial/temporal resolution in some 
MRI applications, parallel MRI (pMRI) acquisition techniques with multiple 
coils acquisition have emerged since the early 1990 's as powerful 3D imag- 
ing methods that allow faster acquisition of reduced Field of View (FOV) 
images. In these techniques, the full FOV image has to be reconstructed 
from the resulting acquired undersampled /c-space data. To this end, sev- 
eral reconstruction techniques have been proposed such as the widely-used 
SENSE method. However, the reconstructed image generally presents arti- 
facts when perturbations occur in both the measured data and the estimated 
coil sensitivity maps. In this paper, we aim at achieving good reconstructed 
image quality when using low magnetic field and high reduction factor. Un- 
der these severe experimental conditions, neither the SENSE method nor 
the Tikhonov regularization in the image domain give convincing results. To 
this aim, we present a novel method for SENSE-based reconstruction which 
proceeds with regularization in the complex wavelet domain. To further en- 
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hance the reconstructed image quality, local convex constraints are added in 
the regularization process. In vivo experiments carried out on Gradient-Echo 
(GRE) anatomical and Echo-Planar Imaging (EPI) functional MRI data at 
1.5 Tesla indicate that the proposed algorithm provides reconstructed images 
with reduced artifacts for high reduction factor. 

Key words: pMRI, SENSE, Regularization, Wavelet Transform, Convex 
Optimization. 



1. Introduction 

Reducing global imaging time is of main interest in neuroimaging, in 
particular in the study of brain dynamics using functional MRI (fMRI). Ac- 
tually, a short acquisition time lets one improve the spatial/temporal res- 
olution of acquired fMRI data, which leads to a more efficient statistical 
analysis. In addition, by reducing the global imaging time, some additional 
artifacts caused by the patient motion can be avoided. For this reason, par- 
allel imaging systems have been developed: multiple receiver surface coils 
with complementary sensitivity profiles located around the underlying ob- 
ject are employed to simultaneously collect in the frequency domain (i.e the 
/c-space), data sampled at a rate R times lower than the Nyquist sampling 
rate along at least one spatial direction, i.e. the phase encoding one. There- 
fore, the total acquisition time is R times shorter than with conventional 
non parallel imaging. A reconstruction step is then performed to build a 
full Field of View (FOV) image by unfolding the undersampled measured 
ones. This reconstruction is a challenging task because of the low Signal to 
Noise Ratio (SNR) in parallel MRI (pMRI) caused by aliasing artifacts re- 
lated to the undersampling rate, those caused by noise during the acquisition 
process and also the presence of errors in the estimation of coil sensitivity 
maps. The Simultaneous Acquisition of Spatial Harmonics (SMASH) [2] was 
the first reconstruction method introduced by Sodickson and Manning in 
1997, operating in the fc-space domain. It uses a linear combination of pre- 
estimated coil sensitivity maps to generate the missing phase encoding steps. 
Some other /c-space based reconstruction techniques have also been proposed 
like GRAPPA (Generalized Autocalibrating Partially Parallel Acquisitions) 
[3], a generalized implementation of variable-density AUTO-SMASH (VD- 
AUTO-SMASH) [4], which improves the AUTO-SMASH [5] method. The 
main difference between SMASH and these methods is that SMASH needs 
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a separate coil sensitivity map estimation, while the other methods need 
the acquisition of few additional /c-space lines to derive coil sensitivity maps 
without a reference scan. This is the reason that motivates the autocali- 
brating qualification of such methods. Indeed, GRAPPA may be preferred 
to non-autocalibrated methods when accurate coil sensitivity maps may be 
difficult to extract, or when reference scans are not possible to obtain at the 
right spatial resolution. For instance, in lung and abdomen imaging, the nu- 
merous inhomogeneous regions with a low spin density make inaccurate the 
estimation of the sensitivity information. However, all the reported methods 
may suffer from phase cancellation problems, low SNR during the acquisition 
process and limited reconstruction quality. 

In [6], an alternative reconstruction method called SENSE (Sensitivity 
Encoding) has been introduced. It is a two-step procedure relying first on 
a reconstruction of reduced FOV images and second on a spatial unfolding 
technique, which amounts to a weighted least squares estimation. This tech- 
nique requires a precise estimation of coil sensitivity maps using a reference 
scan (usually a 2D Gradient-Echo (GRE)). 

To the best of our knowledge, in actual clinical daily routines, only GRAPPA 
and SENSE are offered by scanner manufacturers (or variants like mSENSE 
by Siemens). However, SENSE which is investigated in this work, remains 
the most frequently employed technique. In addition to neuroimaging, many 
other applications like cardiac imaging [7] benefit from the enhanced image 
acquisition capabilities of this method. For a general overview of reconstruc- 
tion methods in pMRI see [8]. 

Actually, SENSE is often supposed to achieve an exact reconstruction 
in the case of non-noisy data and perfect coil sensitivity maps knowledge, 
which is also true for all above mentioned methods. However, in practice, 
noisy data and inaccuracies in the estimation of coil sensitivity maps make 
the reconstruction problem ill-conditioned. To overcome this limitation, some 
alternatives have been proposed like the optimization of the coil geometry [9] 
and the improvement of coil sensitivity profile estimation [10]. But, when low 
magnetic field intensities (up to 1.5 Tesla) are used, high reduction factors 
are considered as unfeasible since the reconstructed images are affected by 
severe aliasing artifacts, and thus, a reduction factor value of two is often 
considered as the highest possible reduction factor at this field intensity. If 
one wants to further improve the spatial resolution to avoid filtering that 
blurs activation, to improve the temporal resolution in fMRI or to reduce 
the global imaging time in clinical MRI, it would be at the expense of the 
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reconstructed image quality and/or potentially the reliability of statistical 
analysis in fMRI and brain activation detection. 

As the image reconstruction problem is ill-posed, several works have dealt 
with regularization methods [1, 11, 12] to better estimate full FOV images, 
most of them operating in the image domain. For an introductory survey of 
linear inverse problems in application to pMRI reconstruction see [13]. 
Standard regularization can improve the reconstructed image quality when 
the experimental conditions are not too severe. However, under a low mag- 
netic field intensity with a high reduction factor, the problem becomes more 
difficult because of the low SNR during the acquisition process and poor es- 
timation of coil sensitivity maps. In fact, in this paper, we are interested in 
a pMRI reconstruction problem under a low magnetic field with high reduc- 
tion factor. We will show later that standard regularization techniques may 
not provide convincing solutions to this problem in such severe experimen- 
tal conditions. Indeed, these methods assume Gaussianity of the solution 
in the spatial domain, which is not a quite good approximation since the 
empirical histogram of an MRI image voxels may be multimodal. So, it is 
mandatory to design more sophisticated methods to reduce aliasing artifacts 
so as to ensure an acceptable reconstructed image quality. In our previous 
work [1] , the appealing properties of wavelet transforms in generating sparse 
representations of regular images have been exploited to build a regularized 
reconstruction technique in the wavelet domain that reduces artifacts. This 
work is extended in this paper to further improve the reconstructed image 
quality, in particular by exploiting the artifact features. We will show that 
the choice of wavelet representations has also been motivated by the ability 
to employ tractable statistical models and the emergence of fast convex non 
differentiable optimization methods. The main contributions of this paper 
are the following: 

• the design of a novel wavelet-based inversion technique allowing us to 
incorporate useful constraints on the images to be reconstructed, 

• the proposition of an efficient convex optimization algorithm that min- 
imizes the related regularized criterion. 

This paper is organized as follows. In Section 2, we give a brief overview of 
the SENSE method and its regularized version in the space domain. Sec- 
tion 3 describes the proposed wavelet-based regularization approach and the 
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associated convex optimization algorithms. Results conducted on anatomi- 
cal GRE and functional EPI data are provided and commented in Section 4. 
Finally, some conclusions and perspectives are drawn in Section 5. 



2. Background 

Although pMRl is a 3D imaging technique, it proceeds with a slice by slice 
acquisition to get all the 3D volume, which simplifies the problem to a two- 
dimensional one. In what follows, we are only interested in reconstructing a 
given slice (2D image). 

2.1. Basic SENSE reconstruction 

An array of L coils is employed to measure the spin density p into the object 
under investigation.^ The signal dg received by each coil £ {1 < i < L) is the 
Fourier transform of the desired 2D field p weighted by the coil sensitivity 
profile Si, evaluated at some locations k = [ky, k^Y in the fc-space. The 
received signal di is therefore defined by the sampling scheme 



where n{k) is a realization of an additive zero-mean Gaussian noise and 
r = [y, xY is the spatial position in the image domain. For the sake of 
simplicity, a Cartesian coordinate system is generally adopted in whole brain 
imaging. 

In its simplest form, SENSE imaging amounts to solving a one-dimensional 
inversion problem due to the separability of the Fourier transform. Note how- 
ever that it extends to a two-dimensional inversion problem for 3D imaging 
sequences like EVI (Echo Volume Imaging) [12] where undersampling occurs 
in two directions of the /c-space. 

Let Ay = ^ be the sampling period where Y is the size of the FOV along the 
phase encoding direction, let y be the position in the image domain along 
the same direction, x the position in the image domain along the frequency 
encoding direction and R < L the reduction factor. A 2D inverse Fourier 
transform allows us to recover the measured signal in the spatial domain. By 
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accounting for the undersampling of the fc-space hj R, (1) can be reexpressed 
in the following matrix form: 



where 



Sir 



A 



d{r) = S{r)p(r) + n(r) 



/ si{y,x) ... si{y + (R- l)Ay,x) \ 
\sLiy,x) ... SLiy+iR-l)Ay,x) J 



iLxR 



(2) 



A 



\^ piy+{R-l)Ay,x) J 
I ni{y,x) \ 



d{r, 



A 



\ dL{y,x) J 



and n(r 



A 



n2{y,x) 



\ nL{y,x) J 



(3) 



In (2), {n{r))r is a sequence of circular zero- mean Gaussian complex- 
valued random vectors. These noise vectors are i.i.d. and spatially inde- 
pendent with covariance matrix [6, 14] of size L x L. In practice, is 
estimated by acquiring data from all coils without radio frequency pulses, 
and its generic entry \l/(£i,£2) corresponding to the correlation between the 
two coils ii and £2 is given by: 



^nEfa.x) se^{y,x)sl{y,x) 

(E(.,x) My,xW){J2(y,^) My,x 



V(£i,£2)e{i,...,L}^ 



(4) 

where (■)* stands for the complex conjugate. 

Therefore, the reconstruction step consists in inverting (2) and recovering 
p(r) from d{r) at spatial positions r = {y,xy. Note that the data {d£)i<i<L 
and the unknown image p are complex- valued, although |p| is only considered 
for visualization purposes. 
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A simple reconstruction method also called the SENSE approach [6], is 
based on the minimization of the Weighted Least Squares (WLS) criterion. 
The objective is to find a vector PwLs(^) each spatial location r such that: 

PwLs(^) =argniin JwlsIpIt-)) 

p(r) 

= (5"(r)*'i5(r))*5^(r)*-M(r) (5) 

where JwLs{p{r)) =\\ d{r) - S{r)p{r) {-Y (resp. {-f) stands for 

the transposed complex conjugate (resp. pseudo-inverse) and, || • = 
^{■Y^-\-) defines a norm on (L . 

In practice, the performance of the SENSE method is limited because of 
the presence of i) distortions in the measurements d{r), ii) possible ill- 
conditioning of S{r) for some locations r and iii) the presence of errors in the 
estimation of S{r) mainly in the center and the brain air interface. These 
undesirable effects are illustrated in Fig. 3 (SENSE reconstruction) which 
shows aliasing artifacts in the reconstructed images for two values of the re- 
duction factor R = 2 and R = A. To increase the stability of the solution 
and reduce the artifacts, a Tikhonov regularization is usually applied. 

2.2. Tikhonov regularization 

As shown in [11, 15], a Tikhonov regularization [16] can be applied to solve 
such an ill-posed inverse problem. The regularization process typically con- 
sists in computing PpwlsI'") minimizer of the following Penalized 
Weighted Least Squares (PWLS) criterion: 

JpwLs{p{r)) = JwLs{p{r)) + k \\ p{r) - p^r) \\]^ . (6) 

where Ir is the i?-dimensional identity matrix. The regularization parameter 
K > ensures a balance between the closeness to the data and the penalty 
term, which controls the deviation from a given reference vector p^{r). The 
solution ppwLs(^) admits the following closed-form expression: 

PPWLS(^) = Pr(^) + 

(5^(r)*-i5(r) + kI S"" {r)^-\d{r) - 5(r)p,(r)). 

(7) 

Note that the accuracy of the solution depends on the reference vector p^{r) 
and the choice of the regularization parameter k. On the one hand, if k 
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is close to zero, we tend to the WLS solution PwLS- ^^e other hand, 
when K takes high values, the penalty term dominates and the solution tends 
to the reference pj.. To select a relevant value for k, some works have pro- 
posed an efficient choice of the regularization parameter [17, 18, 12] by using 
for instance the discrepancy principle or the L-curve technique in order to 
enhance the reconstructed image quality. However, it is worth noticing in 
practice that any quadratic regularization procedure tends to introduce blur- 
ring effects, except for specific values of k. An example of blurring effects 
can be seen in Fig 3 (Tikhonov reconstruction) which displays reconstructed 
full FOV images using Tikhonov regularization. To overcome this limitation, 
one usually resorts to edge-preserving penalty terms [19, 20, 21]. Gener- 
ally, these penalty terms are applied in the image domain and make the 
regularization more efficient by limiting blurring effects and preserving the 
image boundaries. However, the introduction of these terms may lead to 
a non-differentiable optimization problem which is not always easy to solve 
numerically. Here, we propose to apply an edge-preserving regularization by 
proceeding in the Wavelet Transform (WT) domain [22]. We show that it 
is then possible to develop a fast and efficient algorithm to solve the related 
convex non-differentiable optimization problem. 

3. Regularization in the WT domain 

3.1. Motivation 

When carefully analyzing SENSE-based reconstructed images (see Fig. 3 
for R = 4), well spatially-localized artifacts appear as distorted curves with 
either very high or very low intensity. Consequently, we propose to look for 
an image representation where these localized transitions can be easily de- 
tected and hence attenuated. In this respect, the WT has been recognized 
as a powerful tool that enables a good space and frequency localization of 
useful information [22]. In the literature, many wavelet decompositions and 
extensions of these decompositions have been reported offering different fea- 
tures in order to provide sparse image representations. We can mention, for 
example, decompositions onto orthonormal dyadic wavelet bases [23] includ- 
ing the Haar transform [24] as a special simple case, decompositions onto 
biorthogonal dyadic wavelets [25], M-band wavelet representations [26] and 
wavelet packet representations [27]. 

An appealing property of the resulting decomposition is that the sta- 
tistical distribution of the wavelet coefficients can be easily modelled in a 
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realistic way. Hence, the Bayesian framework can be adopted to build an 
efficient reconstruction procedure. Wavelet decompositions have been there- 
fore investigated in a number of works in image processing like denoising 
[28, 29, 30, 31] and deconvolution [31, 32, 33] problems. In medical imaging, 
wavelet decompositions have also been investigated for instance in denoising 
[34, 19, 35], coil sensitivity map estimation [36] and encoding schemes [37, 38] 
in MRI, activation detection in fMRI [39, 40, 41], tissue characterization in 
ultrasound imaging [42] and tomographic reconstruction [43]. 

3.2. Definitions and notations 

In what follows, T stands for the WT operator. It corresponds to a 
discrete decomposition onto a separable 2D M-band wavelet basis performed 
over jmax resolution levels. The objective image p of size Y xX can be viewed 
as an element of the Euclidean space with K = Y x X endowed with 
the standard inner product (■ | ■) and norm || ■ ||. We recall that here we are 
interested in reconstructing one slice (2D image). For this reason, only 2D 
WT operators are investigated. In this context, the following notations are 
introduced. 

Definition 3.1 

Let {ek)i<k<K be the considered discrete wavelet basis of the space . The 
wavelet decomposition operator T is defined as the linear operator: 

T:C'<^C^ (8) 
{{p \ ek))i<k<K- 

The adjoint operator T* serving for reconstruction purposes is then defined 
as the bijective linear operator: 

T* : C-^^ ^ (9) 

K 

{(k)l<k<K ^ CfcCfc- 
k=l 

The resulting wavelet coefficient field of a target image function p is defined 
by C = ((Ca,fc)i<fc<i^,n,ax' (Co,i,fc)i<j<jmax,i<fc<x,) wficrc Kj = K M-^^ IS the 
number of wavelet coefficients in a given subband at resolution j (by assuming 
that Y and X are multiple of M^'^^'^) and the coefficients have been reindexed 
in such a way that C,a,k denotes an approximation coefficient at resolution level 
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jmax and Co,j,k denotes a detail coefficient at resolution level j and orientation 
o e O = {0,...,M - 1}2 \ {(0,0)}. In the dyadic case (M = 2), there 
are three orientations corresponding to the horizontal, vertical or diagonal 
directions. Note that, when an orthonormal wavelet basis is considered, the 
adjoint operator T* reduces to the inverse WT operator and the operator 
norm ||T|| of T is equal to 1. 

3.3. Reconstruction procedure 

An estimate of the objective image p will be generated through the re- 
construction wavelet operator T*. Let ( be the unknown wavelet coefficients 
such that p = T*(. We aim at building an estimate ( of the vector of coeffi- 
cients ( from the observations d. To this end, we derive a Bayesian approach 
relying on suitable priors on the wavelet coefficients. 

3.3.1. Likelihood 

Since the noise has been assumed i.i.d. circular Gaussian with zero-mean 
and between-coil correlation matrix the data likelihood can be expressed 
as: 

gid I TX) = n 9{div) I p(r)) cx exp (- Jl(p)) (10) 

re{i,...,y}x{i,...,x} 

where p = T*( and 

Jdp) = ^WLs(p(r)) 

re{i,...,y}x{i,...,x} 

3.3.2. Prior 

Let us denote by / the prior probability density function (PDF) of the image 
in the wavelet domain. By analyzing the correlation between the real and 
imaginary parts of the wavelet coefficients, very low values of the correlation 
factors have been found for medical pMRI images. Table 1 illustrates the 
values of the correlation coefficient. 

These observations have motivated the choice of independent priors for 
the real and imaginary parts of the wavelet coefficients, the marginal distri- 
butions of which can be separately studied. For the sake of simplicity, we 
also assume that the real (resp. imaginary) parts of the wavelet coefficients 
are i.i.d. in each subband. Their statistical characteristics may however vary 
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Table 1 : Correlation between real and imaginary parts of wavelet coefficients over 3 reso- 
lution levels. 





Approximation 


Detail 


Slice 5 


J = l 


-0.004 


-0.139 


J =2 


-0.031 


-0.140 


J =3 


-0.026 


-0.153 


Slice 9 


J = l 


-0.111 


-0.077 


J = 2 


-0.117 


-0.032 


J =3 


-0.122 


-0.031 



between two distinct subbands. Furthermore, by looking at the empirical 
distributions of the real and imaginary parts of the considered wavelet co- 
efficients, we have noticed that their empirical histograms are well-fitted by 
a Generalized Gauss-Laplace distribution. The histograms present a single 
mode and their shape vary between the Gaussian and Laplacian ones (see 
Fig. 1). The corresponding PDF is: 




Figure 1: Example of empirical histogram of wavelet coefficients. 



At the coarsest resolution level jmax, the distributions of both the real 
and imaginary parts of the approximation coefficients can be modelled by a 
Gaussian distribution since they belong to a low frequency subband. 
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3.3.3. Bayesian inference 



Based on the prior and the hkehhood given hereabove, the MAP estimator 
is computed by minimizing the following criterion: 

C = argmax (ln/(C) + ln^(d | T*()) 

= arg min Jwt(C) (12) 

where Jwt(C) = MTX) + Jp(C) and 



k=l o60 j=l k=l 



with 



(Re(Ca,fc) - /i^^)^ , (Im(C,,) - 



<^oAUk) = all\Re{C,,,)\ + l^|Re(Co,>)r + aj:^|Im(Co,,fc)| + ^|Im(Co,,fc) 



'2 



(«,?:, «J-) e (M;)2, (/3£^/?]:^) e (M;)2, (/x^^/x^-) e and (aRe,ai^) G . 

Hereabove, Re(-) and Im(-) (or ■^'^ and ■^™) stand for the real and imag- 
inary parts, respectively. It is clear that JTwt is convex. However, the op- 
timization procedure cannot rely on conventional convex optimization tech- 
niques like the pseudo-conjugate gradient: although J7l is differentiable with 
a Lipschitz-continuous gradient, jTp is not differentiable. This is a main diffi- 
culty which is frequently encountered in inverse problems involving sparsity 
promoting priors [44, 45]. Therefore, we propose to apply a generalized form 
of the iterative optimization procedure developed in [46] which extends the 
one in [47] and is based on the forward-backward algorithm. 

3.4. Forward-Backward optimization algorithm 

The minimization of Jwt is performed by resorting to the concept of 
proximity operators [48], which was found to be fruitful in a number of 
recent works in convex optimization [46, 49]. In what follows, we give a 
brief overview of this key tool for solving our optimization problem. 
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Definition 3.2 [48] 

Let Tq{x) be the class of lower semicontinuous convex functions from a 
separable real Hilbert space x to ] — oo, +0x3] and let ip e ro(x). For every 
X e x? the function </?(■) + || ■ — x|p/2 achieves its infimum at a unique point 
denoted by prox^x. The operator prox^ : x ~* X '^s the proximity operator 

Ofifi. 

Here is an example of proximity operators wliicli will be used in our approach. 
Example 3.1 

Consider the following function: 

(/?:M^R (14) 



with a G IR+, P G M.^ and /i G M. The associated proximity operator is given 
by: 

P + 1 

where the sign function is defined as follows: 



prox^^ = ^-y-^ max||^-/i| -a,0} + /i (15) 



Ve G M, sign(0 



1 if^>0 
-1 otherwise. 



In this work, as the observed data are complex-valued, we generalize the 
definition of proximity operators to a class of convex functions defined for 
complex-valued variables. 
For the function 

$: ^] - oo,+oo] (16) 
X ^ 0^^(Re(x)) + 0^°'(lm(x)), 

where 0^*^ and 0^™ are functions in ro(M^) and Re(x) (resp. Im(x)) is the 
vector of the real parts (resp. imaginary parts) of the components of x G C^, 
the proximity operator is defined as: 

prox^:C^^C^ (17) 
X t-^ prox^Rc(Re(a;)) -|- zprox^im(Im(x)). 

Here is an example of proximity operator for a function of a complex-valued 
variable. 
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Example 3.2 

Consider the following function: 



$ : C^M (18) 

y 



^Re 



olm 

ai-|lm(e-/i)| + ^|Im(e-/i)P 



with e (M+)2, (/^^^/^i"^) G (M;)^ and G C. The associated 

proximity operator is: 



prox^e = "''^"j,'^!'! 1 ^" max{|ReK - m)| - a■'^ 0} 



+ ' "^yfr" - f )l - 0} + ^'■ (19) 

Based on these definitions, and by extending the algorithm in [46] to the 
complex case, a minimizer of JTwt can then be iteratively computed according 
to Algorithm 1. Note that in this algorithm, the expressions of prox^^^^ and 
prox^^ ij,^^ at each iteration n are provided by Example 3.2. It can also be 
noticed that A„ and 7„ respectively correspond to relaxation and step-size 
parameters. 

For 3D volume reconstruction we iterate over slices, and for functional 
data, we then iterate over volumes independently (see discussion in Sec- 
tion 5). 

3.5. Convergence of Algorithm 1 

For every r G {1, . . . , Y} x {1, . . . , X}, let 6'r > be the maximum eigen- 
value of the Hermitian positive semi-definite matrix S'^(r)^~^S'(r) and let 
9 = maxre{i,...,y}x{i,...,js:} > 0. To guarantee the convergence of Algorithm 
1, the step-size and relaxation parameters are subject to the following con- 
ditions: 

Assumption 3.3 

(i) inf„>o7„ > and sup„>o7n < ^rp, 

(ii) inf„>o Xn > and sup„>o -^n < 1- 
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Algorithm 1 2D-slice wavelet-based regularized reconstruction 

Let (7n)n>o and {Xn)n>o be sequences of positive reals. 

1: Initialize C^^^- Set n = 1, e G (0, 1) and J^^^ = 0. 
2: repeat 

3: Reconstruct the image by setting p*^") = T*(^'^\ 
4: Calculate the image such that: 

VrG{l,...,F}x{l,...,X}, 

*i(")(r) = 25^(r)*-i (5(r)pW(r) - d(r)), 

where the vector u^"'\v) is defined from m*^") in the same way as p(r) 

is defined from p (see (3)). 
5: Determine the wavelet coefficients v^""^ = Tu^"'^ of 
6: Update the approximation coefficients of the reconstructed image: 

7: Update the detail coefficients of the reconstructed image: 
VogO,Vj e {l,...,j^ax},VA;GK,-, 

/■("+!) _ An) , X / //-(n) _ (n) n _ >(n) \ 

'^o,j,k ~ So,j,fc \^P^"-^7„$oj I'soj.fc in'Jo,j,k) ^o,j,kJ- 

8: Compute = Jwt(C^"^)- 

9: n ^ n + 1 
10: until - < 

11: return = T*C^"^ 
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More precisely, the following result can be shown: 

Proposition 3.4 Under Assumption 3.3, the sequence (C*-"''')n>o built when 
iterating Steps 3 to 9 of Algorithm 1 converges linearly to the unique solution 
C to Problem (12). 

Proof: See Appendix A. 

Results obtained with the proposed iterative method are provided and 
discussed in Section 4. In these results, it can be noticed that some artifacts 
remain in the reconstructed full FOV images (see the image obtained with 
Algorithm 1 in Fig. 4) because they present extremely large volumes. In 
the next section, we therefore present an extension based on additional con- 
straints that enable a very accurate reconstruction and cancellation of such 
artifacts. 

3.6. Constrained wavelet- based regularization 

We propose to extend our approach by incorporating an additional con- 
straint in the method described hereabove in order to better regularize arti- 
fact regions. 

3.6.1. New optimality criterion 

As artifacts appear as curves with very high or very low intensity, we 
propose to set local lower and upper bounds on the image intensity values in 
artifact areas. These bounds define the nonempty closed convex set: 

C = {p G I Vr G {1, . . . , r} X {1, . . . , X}, p(r) G a} (20) 

where the constraint introduced on the range values at position r G {1, . . . , Y} x 
{1, . . . , X} is modelled by 

a = {e e c I Re(o G im(o g fn, (21) 

with I^'^ = [Im'in,r5 Imax.r] ^^^^ I^™ = [l[jjin,r ' Imax,?-] • 

When taking into account the additional constraint, the constrained cri- 
terion in (12) becomes: 

Jcwt(C) = MO + Jp(C) + ^c.(C), (22) 

where 

C* = {C G C^' I T*C G C} 
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and ic* is the indicator function of the closed convex set C* defined by 



if C G C* 
-oo otherwise. 



The minimization problem can then be rewritten as: 

Find mmMO + JpiO+tciO- (23) 

3.6.2. Computation of the new estimator 

Conceptually, in order to solve the minimization problem in (23), the 
forward-backward iteration in (36) (in Appendix A) has to be updated ac- 
cording to: 

^(n+i) ^ ^(n) ^ (^prox^^^^^^^^ W _ 7„V:ri(C("))) - C^"^) • (24) 

The main difficulty here is that the proximity operator of 7.„j7p + zc* does not 
have a closed form. However, from the definition of the proximity operator, 
we get 

VC e C^, prox,„^^+,^, (C) = arg min + J' dO (25) 

where 

^'c(-) = ^l|--Cir + ^c.(-)- (26) 

Although prox^^ does not take a simple expression, the proximity oper- 

ator of 7nJ7p is expressed by (35) (in Appendix A) and the proximity operator 
of J' is easily determined. Indeed, from simple calculations, we have 



VC'GC^, prox^,,(0=i^c. 



(27) 



where Pc* is the projection onto the convex set C* . In turn, provided that the 
considered wavelet basis is orthonormal, the projection onto C* of C,' G 
is obtained by performing the wavelet decomposition of the projection of 
p' = T*(' onto C. The latter projection is 

Pcip') = {PcAp'ir))^^^^_^^^^^_^^ (28) 
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where, for every r G {1, . . . , Y} x {1, . . . , X}, 

Cax,. ifRe(0>eax,. (29) 
^ otherwise, 

a similar expression being used to calculate Im(Pc'^(^)). 

Knowing prox^^j.^ and pioxj,^, pTox^^jp_^^^X can be computed in an it- 
erative manner by solving the optimization problem in (25) with the Douglas- 
Rachford algorithm [50, 51]. More precisely, we have: 

Proposition 3.5 

Set r]^^^ E C^"" and construct for all m E N: 
j^(m+|) _ proxj-z^r/^™^ 

^(m+i) ^ ^(m) ^ r(prox^^^^(2r/(™+5) - r/M) - r/('"+i)), 

where r g]0,2[. Then, {ri^"^^^^)rn.eN converges to prox^^j-^_,_j^^C- 

Inserting this extra iterative step in the forward-backward algorithm and 
using the expressions of prox^^j^^ and pioxj,^ lead to Algorithm 2. 
At iteration n, Mn is the number of times the Douglas-Rachford step is run. 
According to the convergence analysis conducted in [50, Prop. 4.2], if M„ 
is chosen large enough and Assumption 3.3 holds, iterating Steps 3 to 16 of 
the algorithm guarantees the convergence to the unique solution to Problem 
(23). Note however that [50, Prop. 4.2] does not provide a practical guideline 
for the choice of M„. The practical rule we chose is explained in Section 4.1. 
The improvements resulting from the use of this constrained approach are 
illustrated in the next section. 
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Algorithm 2 Constrained 2D-slice wavelet-based regularized reconstruction 
Let (7n)n>o and (A„)„>o be sequences of positive reals, let 
{Mn)n>o be a sequence of positive integers and set r g]0, 2]. 

1: Initialize C^^^- Set n = 1, e G (0, 1) and J^^^ = 0. 
2: repeat 

3: Same as for Algorithm 1 
4: Same as for Algorithm 1 
5: Same as for Algorithm 1 

6: Initialize the Douglas-Rachford algorithm by setting 

^(n,0) _ ^(n) _ ^^^(")_ 

7: Douglas-Rachford iterations: 

• for m = to Mn — 1 do 

9: Compute ry("'"*+2) = p^. (^1 ^^~) ' 

10: Update the approximation components of r/('^''"): 

V/C t Ji^jmax; '7a,fc — Va,k +^ lPr0^7n*a l^'^a.fc Va,k ) 

'la,k ) ' 

11: Update the detail components of r^^"''"): 
VoeO,Vj e {l,...,j„,ax},VfceK,-, 

12: If r/("'™+i) = goto 14. 

• end for 

14: Update the wavelet coefficients of the reconstructed image: 

^(n+l) ^ ^(n) _^ A„(r/("'"*+5) — 

15: Compute = Jcwt(C^"^)- 
16: n ^ n + l. 
17: until - < 

18: return p("' = r*C^"^ 
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4. Experimental results 

Experiments have been conducted on real data sets containing 256 x 
256 X 14 Gradient-Echo anatomical and 64 x 64 x 30 EPI functional images 
with respectively 0.93 x 0.93 x 8 (mm^) and 3.75 x 3.75 x 3 (mm^) spatial 
resolution. For Gradient-Echo data, we have TE/TR = 10/500 (ms) with 
BW=31.25 (kHz), while for EPI data we have TE/TR = 30/2400 (ms) and 
BW=125 (kHz). These images have been acquired using reduction factors 
R = 2 and i? = 4 to show the effectiveness of our approach on image re- 
construction after a significant reduction of the global imaging time. The 
employed scanner was Signa 1.5 Tesla GE Healthcare with an eight-channel 
head coil. The standard deviation (j„ of the acquisition noise was estimated 
based on the background region of the measured images where there is no 
signal corresponding to the region of interest (in our case the head). In order 
to proceed to an objective assessment of the reconstruction accuracy, a full 
FOV image prcf is available through a non accelerated acquisition for i? = 1. 

4.1. Anatomical data 

In this experiment, dyadic (M = 2) Symmlet orthonormal wavelet bases 
[52] associated with filters of length 8 have been used over j^ax = 3 resolution 
levels. For the wavelet coefficients, the priors described in Section 3.3.3 
have been employed. The related hyper-parameters have been estimated 
based on a reference image pref using the maximum likelihood estimator. A 
couple of hyper-parameters (a,/9) is estimated for real/imaginary parts of 
each subband. Algorithm 1 was then used to reconstruct full FOV images. 
For the sake of simplicity, constant values along the algorithm iterations have 
been adopted for the relaxation parameter A„ and the step-size 7„. It has 
been noticed experimentally that = 1 is the optimal value of the relaxation 
parameter in terms of convergence rate (see Fig. 2). 

A value of the step-size closed to the allowed maximum value in Assump- 
tion 3.3 has also been observed to be beneficial to the convergence rate. 
After computing the constant 6 related to the considered sensivity map, 7„ 
was thus chosen equal to 1.99/^ = 12.83. 

The algorithm has been stopped when the criterion no longer signifi- 
cantly varies, by choosing e = 10^^ in Algorithm 1. For different values 
of A, Fig. 2 illustrates the evolution of the optimization criterion w.r.t the 
iteration number when reconstructing a 2D-slice. Comparison of the curves 
shows that A^ = 1 gives the fastest convergence. One can also consider that. 
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Figure 2: Convergence speed of the optimization algorithm w.r.t. the choice of the relax- 
ation parameter A for jmax — 3. 



after about 20 iterations, which correspond to 160 seconds of execution time 
using Matlab 7.7 on an Intel Core 2 (3 GHz) architecture, the minimum was 
reached with A„ = 1, so providing the optimal MAP solution. 
Note that accelerated algorithms such as TWIST or FISTA have been re- 
cently proposed in the literature [53, 54] for minimizing the same optimality 
criterion. However, for the considered application, it has been observed no 
significant improvement of the convergence profile by using these methods. 
Fig. 4 shows reconstructed full FOV anatomical images using the proposed 
approach (Algorithm 1) for R = 2 and R = A. To compare this reconstruc- 
tion with the Tikhonov one, we can refer to Fig. 3. Note that in Tikhonov 
regularization, the regularization image pr was chosen as a mean image based 
on the basic-SENSE reconstruction, which contains the mean value of the sig- 
nal of interest within the brain region. The regularization parameter k, was 
manually fixed. A comparison with the basic-SENSE reconstruction can also 
be made from Fig. 3. 

It can be observed that the aliasing artifacts in the basic-SENSE recon- 
structed image are smoothed. They are completely removed in the areas 
where they were not actually very strong. When observing the reconstructed 
image using Tikhonov regularization, we see that it suffers from blurring ef- 
fects (Fig. 3). These drawbacks do no longer exist in the WT regularized 
image in Fig. 4 where a good reconstruction is performed in the brain area. 

If we compare the proposed approach with basic-SENSE and Tikhonov 
reconstructions from a quantitative viewpoint, clear improvements in the re- 
construction process are also obtained. The quantitative measure we used is 
the Signal-to- Noise Ratio (SNR). Using the available reference image pref and 
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R = 2 



R = A 



Algorithm 1 



Algorithm 2 




Figure 4: Two reconstructed slices using Algorithm 1 and Algorithm 2 for R — 2 and 
R = A. 
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for the reconstructed image p obtained with any reconstruction approach, it is 
possible to evaluate the signal-to-noise ratio: SNR = 201og]^Q (||pref ||/||Prof — 
pII). Table 2 gives the SNR values corresponding to the basic-SENSE, 
Tikhonov regularization and the proposed WT regularization for different 
slices of the anatomical brain volume. 



Table 2: SNR evaluation for reeonstructed images using different methods. 





SNR (dB) 




SENSE 


Tikhonov regu- 


WT regulariza- 


Constrained WT 






larization 


tion 


regularization 


Slice 1 


14.34 


14.48 


14.67 


15.53 


Shce 2 


11.53 


11.73 


13.72 


15.13 


Slice 3 


12.96 


13.37 


14.02 


14.12 


Slice 4 


9.22 


9.58 


10.16 


13.10 


Slice 5 


11.49 


11.88 


12.06 


12.25 


Slice 6 


9.67 


9.80 


10.11 


11.22 


Slice 7 


11.04 


11.26 


11.52 


12.00 


Shce 8 


12.19 


12.60 


12.92 


13.60 


Shce 9 


13.74 


13.28 


14.29 


15.66 



Note that reconstructed images here correspond to slices of the middle 
region of the brain volume for which full FOV images contain a large area of 
signal of interest w.r.t the image background. 

To improve reconstruction, the constrained algorithm presented in Sec- 
tion 3.6 was applied with the Symmlet 8 wavelet basis. The parameters 
A„ and 7„ have been set to the same values as for the unconstrained case, 
while r was fixed to 2 as it was practically observed that this value gives 
the best convergence rate for the underlying Douglas-Rachford iterations. 
In practice, the value of M„ was defined as the minimal integer value such 
that 1 ^^ ' "r^^flT'^^?^ — -| < 10~^, which results in about 4 iterations of the 

I rq{n,l\in—Z) I 1 

Douglas-Rachford algorithm. A morphological gradient [55] was used to de- 
tect artifact regions on which we apply the additional convex constraint. The 
upper and lowed bounds which define the convex sets Cr (Eq. (21)) were 
computed based on a morphological opening and closing applied to the basic- 
SENSE reconstructed image in order to discard very low and very high in- 
tensities. Reconstructed anatomical images using this constrained approach 
is displayed in Fig. 4. 
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Usual inspection of the results shows that the surviving artifacts in Fig. 3 
have now been removed. A comparison with basic-SENSE, Tikhonov and 
unconstrained WT regularizations can be made from Figs. 3 and 4. Notice- 
able improvements in terms of SNR are also obtained that can reach 3.6 dB 
w.r.t basic-SENSE, 3.4 dB w.r.t Tikhonov regularization and 2.94 dB w.r.t 
unconstrained WT regularization (see Table 2). Finally, it is important to 
note that the performance of our constrained reconstruction also depends on 
the efficiency of the prior detection (based on morphological analysis) of the 
distorted area on which the convex constraints are applied. 

4-1.1. Choice of the wavelet basis 

In this section, we study how the choice of the wavelet basis may influence 
the reconstruction performance. For comparison purposes, we present the 
results obtained with four different wavelet bases: dyadic Symmlet 8, dyadic 
Daubechies 8, dyadic Haar and Meyer with M = 4 bands [56]. In Fig. 5, 
reconstructed images using the different wavelet bases with jmax = 3 are 
displayed. 

Some boundary effects are introduced in the reconstructed images, but 
they are not very strong and do not affect the brain area of interest when 
using the Symmlets, Daubechies and Haar wavelet basis. However, these 
additional artifacts become very important when using the Meyer 4-band 
wavelet basis because of its large spatial support. Hence they drastically 
decrease the SNR of the reconstructed full FOV image. Note also that when 
using the Haar wavelet basis, we introduce some blocking effects caused by 
the Haar wavelet discontinuities that do not occur with the Symmlets and 
the Daubechies bases. Among the latters, Symmlets 8 gives slightly better 
regularization results than Daubechies 8. 

More extensive comparisons using other bases were conducted, and we can 
conclude through this study that Symmlets 8 is the orthonormal wavelet 
basis among the tested ones which achieves the best WT regularization per- 
formance in terms of both SNR and visual quality. 

4.1.2. Choice of the maximum resolution level 

In this section, we focus on the effect of the choice of the maximum resolution 
level jmax in terms of reconstruction quality. 

The impact on reconstructed full FOV images can be emphasized through 
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Symmlcts 8 {SNR = 15.66 dB) 



Daubcchics 8 (SNR = 15.60 dB) 




Figure 5: Reconstructed images using different wavelet bases. 

the difference between reconstructed images using 1, 2 and 3 resolution levels. 
Fig. 6 illustrates the difference between reconstructed images using the first 
and the second (left), the second and the third (right) resolution levels. 

(Image with Jmax — 2) — (Image with Jniax — 1) (Image witli Jmax — 3) — (Image with Jmax — 2) 




Figure 6: Difference between reconstructed images using different resolution levels. 

It appears that the difference between reconstructed images at different 
resolution levels is quite important since it can reach a value of 40 within 
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an intensity range [0,255]. Moreover, this difference may be important in 
distorted areas. Hence, the higher the maximum resolution level, the better 
regularized the artifacts are. However, only slight improvements are obtained 
beyond 3 resolution levels. 

Note that by increasing the maximum resolution level, boundary effects are 
more visible, but they do not affect the brain. 

A maximum resolution level jmax = 3 may therefore be a suitable choice for 
achieving an acceptable reconstructed image. 

4-2. Functional data 

This experiment on functional EPI data of size 64 x 64 was conducted 
using the same wavelet basis and priors. Algorithm parameters (i.e. relax- 
ation and step-size parameters) have been adjusted according to the same 
rules as for anatomical data. The relaxation parameter was set to 1 and 
the step-size parameter was set to 1.99/6' = 20.63 after computation of the 
6 constant. Note that these EPI images have been acquired with no ex- 
perimental paradigm at resting state (eyes closed and subject lying on the 
MRI bed). Fig. 7 illustrates reconstructed full FOV slices using SENSE, 
Tikhonov, Algorithm 1 and Algorithm 2. 

It can be shown that many defective pixels were corrected when using 
the proposed WT regularization in Algorithm 1 without introducing addi- 
tional artifacts, in contrast with Tikhonov regularization. When analyz- 
ing reconstructed images using Algorithm 2, it can be noticed that pixels 
with very high intensity were smoothed, when compared with Tikhonov and 
basic-SENSE reconstructed images. Residual defective pixels belonging to 
distorted areas in SENSE reconstructed image have also been removed due 
to the convex constraint introduced in these areas. Note that the same ap- 
proach as that applied in anatomical experiments has been adopted to detect 
artifact areas and to compute the upper and lower bounds defining the convex 
set • 

4-2.1. Analysis of resting state BOLD signal 

To validate our reconstruction approach, an analysis was carried out on the 
occipital cortex area, where a statistical study is possible even if we have 
resting-state images with no experimental paradigm. Fig. 8 shows the vari- 
ation of the BOLD signal w.r.t. time for images reconstructed with R = 4 
using SENSE, Tikhonov and the proposed approach. The variation of the 
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SENSE 



Tikhonov 



Algorithm 1 



Algorithm 2 




Figure 7: Two reconstructed slices using SENSE, Tikhonov regularization, Algorithm 1 
and Algorithm 2 for i? = 4. 
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BOLD signal of images reconstructed using SENSE with i? = 2 is also given 
and considered as a reference signal since with R = 2, SENSE provides a 
quite reasonable reconstruction quality. 




Figure 8: Time evolution of the BOLD signal in the occipital cortex. 



The four plots indicate that, by using our WT approach (constrained 
wavelet-based regularization), the BOLD signal is smoother than when using 
SENSE or Tikhonov regularization for EPI images acquired with R = 4. In 
fact, the occipital cortex is a brain region which is always activated even 
if there is no experimental paradigm, which justifies the smoothness of the 
BOLD signal in this region. 

5. Discussion and Conclusions 

In this paper, we proposed and tested a novel approach for SENSE re- 
construction based on a regularization in a 2D orthonormal wavelet basis. 
This method reduces aliasing artifacts related to the noise in the acquired 
data, which become critical when using high reduction factors with low mag- 
netic field intensity. Our results on 1.5 Tesla data show improvements in 
reconstructed images compared with basic SENSE and other standard reg- 
ularization methods for both anatomical and functional images. Different 
choices for the wavelet basis and the maximum resolution level have been 
discussed to evaluate the impact of these configurations on our algorithm 
performance. Another potential advantage of the proposed method is that 
it can be implemented with parallel CPU architectures as we can deal with 
each wavelet coefficient subband separately in some steps of the algorithm. 
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Parallelism can also be exploited to reconstruct different slices or, in fMRI, 
different volumes independently. Since our approach is more time-consuming 
than SENSE when reconstructing a 2D-slice, this parallelism let us catch up 
computational time when reconstructing whole brain volumes or temporal 
series of brain volumes. 

A constrained version of our algorithm was also presented in order to achieve 
more accurate reconstructed images. Our results show improvements both 
in visual image quality and in quantitative measures. 

Concerning our future work, a more sophisticated detection of artifact regions 
could improve the obtained results. The incorporation of additional convex 
constraints on artifact regions could also be of great interest. An extension to 
3D wavelet decompositions may also be envisaged to reconstruct whole brain 
volumes at once and probably exploit between-slices spatial dependencies. 
Experimentations of our approach on 3 Tesla data and functional data with 
experimental paradigm may also be important. 



From (13)-(14), it can be seen that Jp is a convex function such that 



This means that Jp is a strongly convex function with modulus 'dh'^ ■ Since 
J7l is a finite convex function, JTwt also is strongly convex. It is thus strictly 
convex and coercive (i.e. lim|j^||^+oo J7wt(C) = +00) and, from standard 
result in convex analysis [57], it can be deduced that Jwt has a unique 
minimizer C,. 



function /: x ^\ ^ cxD; +00] i where x is a Hilbert space, is said strongly convex on 
X with modulus > if there exists some g & To{x) such that f ~ g + ^^^^^ ■ 



APPENDIX A 
Proof of Proposition 3.4 



where 




^1 



min{(2a^j,) \ (2crij^) \ (/3j>)oGO,i<j<j,nax5 (/^j>)oeo,i<j<j,nax}- 
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In addition, J7l is a differentiable function and we have 

VCGC^, V:r.HC) = ^^ + ^^^ = TVJ.(T*C), (31) 

where Jlt{C) = JhiTX). Set p = T*C and u = VJl(p)- It can be then 
deduced from (11) that 

Vr e {1, . . . , F} X {1, . . . , X}, n(r) = 25"(r)*-i (5(r)p(r) - d(r)) . 

where the vector u{r) is defined from u in the same way as p(r) is defined 
from p in (3). Furthermore, for every (' G C'^, 

WVJLTiO-'^JLTiOW < \\T\\\\u-u'\\ (32) 
where u' = VJl{p') and p' = T*('. We have then 
— m'||^= ||tt(r) — u'(r)|p 

re{l,.--,'^}x{l,...,X} 

= 4 Yl ||5^(r)*-^5(r)(p(r)-p'(r))||2 

re{l,...,Y}x{l,...,X} 

<4 ^ ^r'llP(r)-p(r)f 

re{i,...,r}x{i,...,x} 

<4^'IIP-P'f 

<Ae^\\T\\^\\C-Cf- (33) 
Altogether, (32) and (33) yield 

||VJlt(C) - VJlt(C')II < 29\\Tf\\( - C\\ (34) 

which shows that J'lt has a Lipschitz continuous gradient with constant 
2^||T||2. 

Based on these observations and the fact that, 

(Co,,,.) max : 

Prox7nJpC = ((prox^„<i.Xa,fc)i<fe<i^.„.ax ' (prox^„*<,,,Co,i,fc)i<i<j,„ax,i<fe<i^,) , (35) 

the sequence (C'^"''')n>o built by Algorithm 1 can be rewritten under the more 
classical forward-backward iterative form [58]: 

^(n+i) ^ ^(n) + (prox^^^^ (C^") - 7nV Jz.t(C^"^)) - C^"^) (36) 
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and, due to the Lipschitz differentiability of Jlt-, the convergence of the 
algorithm is secured under Assumption 3.3 (see [46, 49]). Furthermore, since 
J7p is strongly convex with modulus di, we have (see [50] and references 
therein) 

Vn > 0, IIC^") - Cll < (l - - Cll (37) 

where 7 = inf„>o7n and A = inf„>o A^. This proves that (C''"'')n>o converges 
linearly to C,. 
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